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FOREWORD 


This study was initiated under the direction of Dr. R.R. Jayroe, Flight Data 
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Marshall Space Flight Center. This work was accomplished via the in-house efforts of Dr. 
R.R. Jayroe and Mr. C.W. Campbell of the same office, and by the efforts of 
Dr. J.F. Andrus of Northrop Services, Inc., Huntsville, Alabama, under contract 
NAS8-21810. Dr. Andrus must be credited for conceiving the original idea used in the 
binary correlation routine, and has continued to contribute to this study since joining the 
faculty of Louisiana State University, New Orleans, in August 1973. 



DIGITAL IMAGE REGISTRATION METHODS BASED UPON 
BINARY BOUNDARY MAPS 

I. INTRODUCTION 


The ability to register or match the ground scene of images in the form of digital 
data serves an important role in the analysis of remotely sensed, multispectral, earth 
observation data and change detection. In the case of multispectral camera data, the 
images from the different camera stations must be registered before a proper analysis can 
be performed. Registration is also necessary when the data are acquired from several 
different sensors. In the case of change detection, where data are acquired from the same 
ground scene but at different time intervals, registration has an important role in 
determining what changes in shape have taken place in the ground scene, and registration 
permits an analysis to determine how the multispectral signatures of various ground 
features change as a function of time. For meteorological applications, registration 
methods based on meteorological satellite imagery can be used to estimate cloud 
velocities for global weather information. 

Because of the typically large volumes of data involved, the use of binary maps in 
registration is recommended for the reasons discussed below. Converting the raw data to 
a binary boundary map typically represents a data compression of 60 to 70 percent of 
the original data. An additional significant compression of data is realized by working 
with sequences of, rather than individual, boundary points, which requires only the 
knowledge of the start scan and column and the length of the sequence. 

If A,(x,y) and A 2 (x+£, y+f) represent the amplitudes of the data in images 1 and 
2 respectively, at scan coordinates x and x+£ respectively, and at column coordinates y 
and y+f respectively, then the average product A, (x,y) A 2 (x+£, y+J) is computed over x 
and y for various combinations of £ and f to determine the scan and column shifts 
necessary to register the two digital images. If this average is applied to the raw data, it 
may also be necessary to remove mean values and normalize in order to produce a 
correlation coefficient with a well defined peak. The use of the raw data appears to have 
some drawbacks as compared to using binary border maps. First, an ambiguity can result, 
since a large negative correlation peak can occur and represent a good match as well as a 
large positive peak. This is in most part due to the variation of the data acquired under 
different environmental conditions; i.e., different sun angles, seasons, atmospheric 

conditions, etc. On the other hand, binary boundary maps are produced from relative 

changes in the data, and the boundaries indicated in the ground scene tend not to change 
with environmental conditions. Since the binary boundary maps contain data that consist 
only of 0’s and l’s, the average product A, (x,y) A 2 (x+£, y+f) will always be positive and 
tends to be sharp so that removal of means and normalization is not necessary. In 

addition, it is possible to compute the average product A, (x,y) A 2 (x+£, y+f) without 

multiplying when binary data are used. Addition can be used to replace multiplication, 
and this reduces computer computation time significantly. 



Section II is a discussion of the production of binary boundary maps, while 
Sections III and IV are concerned with the correlation and registration schemes, 
respectively. Section V contains the results, test cases, and summary 


II. BINARY BOUNDARY MAPS 


The purpose of the binary boundary map is to categorize the digital data 
representing the ground scene into homogeneous areas and boundaries, and an example of 
a boundary map from KRTS digital imagery with some feature identification is shown in 
Figure I . 


The computer program used to generate the boundary map f pically runs about 
13 minutes (IBM-7094 time) on an area of 1000 scans by 255 columns, which is 

approximately 327 samples/second. This 
time includes the amount of time required 
to read and compile the program in addition 
to the calculation time. At present, no 
concerted effort has been made to optimize 
the running time of this program or to 
specialize its output for the registration 
program, but the efforts are in a holding 
i n status for future consideration. 


POSSUM HOLLOW 


TENNESSEE RIVER , 


HOBBS ISLAND 


• HUNTSVILLE MOUNTAIN 


, GREEN MOUNTAIN 


WALLACE MOUNTAIN 


The logic of the program is to work 
with two scans of raw digital data 
simultaneously. Let jAjj represent the 

algebraic value of the data in the kth 
channel of data at scan i and column j as 
shown in Figure 2. If there are n channels of 
data, then k ranges from 1 to n. For each 
location i, j, the average vector component 
distances, which represent various degrees of 
spectroscopic changes within the ground 
scene, are computed using the formulas 
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Figure I. Boundary map from a 
portion of FRTS imagery. 
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Figure 2. Computation of S x and Sy for channel K. 


and are stored in a joint histogram, P(S X , Sy). Dividing by n in equation (1) assures that 
S x and Sy have the same typical range regardless of the number of data channels. Based 

upon previous experience with various types of data, it was found that the size of the 
joint histogram array could be limited in size to a 50 by 50 array such that any data 
element with a value of S x or Sy larger than 50 could be automatically assigned as a 

boundary element. As the raw data are read into the program and the number of 
occurrences of S x and Sy are accumulated in the joint histogram, a new tape is created 

for later use. This tape contains the numbers 1 through 50 2 = 2500 and instead of 
writing S x and Sy on this tape for each data element, one number is written on the tape 

that gives the location of S x and Sy in the joint histogram. The location, 1, of S x and Sy 
in the joint histogram is computed using 


I = (5 1 ) Sjj + Sy 


( 2 ) 


and is a unique one-dimensional representation of the two dimensions S x and Sy. The 

new tape containing a number 1 for each data element eliminates the necessity for 
recalculating S x and Sy for each data element a second time. After all the data elements 

from the data set have been exhausted and the joint histogram is complete, a decision is 
made as to which combinations of S x and Sy are to be considered as boundaries. The 

decision curve for a boundary element is based upon the formula 
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where S x and Sy are the values of S x and Sy -at the mode of the joint histogram and 

ASCN, ACOL, IPOW, and BLIM are input parameters to the program. The input 
parameters allow for a wide variety of decision curve shapes and positions. Nominally 
IPOW = 2 and BLIM = 1 , whereas ASCN is usually equal to ACOL and must be estimated 
based upon experience. The possible values of I or correspondingly S x and Sy are then 

inserted into equation (2) and the variable N(I) in the computer program is set equal to 0 
if equation (2) is not satisfied and is set equal to 1 if equation (2) is satisfied. The tape 
containing I for each d”ta element is then read into the program, and the value of N(I) is 
written out as another tape which is called the binary boundary map and contains only 
the numbers 0 and 1. The use of the intermediate tape containing the I values permits 
the boundary decision curve to be varied without recalculating S x and S y for every data 
element. 


III. CORRELATION SCHEME 


The development work for the correlation scheme was initiated because of a lack 
of a readily available registration program and after a review and discussion of some of 
the most recently published papers on registration, such as References 1 through 3. A 
computer listing of the correlation scheme is provided in Appendix A. The central idea is 
to compute only those correlation delays which are affected by overlapping boundary 
elements for each shift of the window against the picture. 

Before describing the correlation scheme, it is necessary to define the computer 
program variables and terminology for the computer routine. 

Picture A section of a boundary map to be used as the reference data. 

Window A section of a boundary map which is to be registered with the 

picture. The window has one-half as many scans and columns of 
data as the picture. 

NROW The number of data scans in the window. 

NCOL The number of data samples per scan or columns of data in the 

window. Nominally, NROW is equal to NCOL. 
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NCOR (I,J) 
(NROWXNCOL) 

KMAX 

LMAX1 

LMAX 

NW(J) 

K 

L 

IW(K) 

JW(K) 

NSQW(K) 

IP(L) 

JP(L) 

NSQP(L) 

I 


The average product A,(x,y) Aj(x+I, y+J). The value of 
NCOR(I,J) is the number of coinciding boundary points that 
occur when the upper left comer element of the window is 
placed upon the 1th row and Jth column of the picture. 
(I* 1,2, NROW + 1 and J = 1,2 NCOL+ 1.) 

Total number of boundary sequences in the window. 

Total number of boundary sequences in the picture for 
1-1,2, ..., NROW. 

Total number of boundary sequences in the picture for all 1, 
1-1,2, .... (2) NROW. 

Dummy variable for reading binary boundary map data. 
(J * 1,..., NROW for window; J= 1, .... (2)NROW for picture.) 

Number of the Kth window boundary sequence 
(K = 1,..., KMAX). 

Number of the Lth picture boundary sequence 
(L = 1,..., LMAX). 

The value of 1W(K) is the scan number minus one of the Kth 
boundary sequence in the window. 

The value of JW(K) is the column number minus one of the 
beginning of the Kth boundary sequence in the window. 

The value of NSQW(K) is the number of boundary elements 
minus one in the Kth boundary olumn sequence in the 
window. 

The value of IP(L) is the scan number of the Lth boundary 
sequence in the picture. 

The value of JP(L) is the column number of the beginning of 
the Lth boundary sequence in the picture. 

The value of NSQP(L) is the number of boundary elements 
minus one in the Lth boundary column sequence in the picture. 

The value of I is the scan shift in NCOR(U). 
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JL The value of JL is the lower value of the column shift J in 

NCOR(I,J). 

JU The value of JU is the upper value of the column shift J in 

NCOR(l,J). (J = JL, JL + 1,...,JU.) 

NAD The value of NAD is the number of coinciding boundary points 

at a given I and J for the Kth window boundary sequence and 
the Lth picture boundary sequence. 

JM The value of JM controls the value of NAD*. JM causes NAD to 

decrease as the Kth and Lth boundary sequences decreasingly 
overlap because of an increase in J from JL to JU. 

NU For a given I, the value of NU is the column number of the end 

of the Lth boundary sequence in the picture. 

NSPAN For a given I, the value of NSPAN is the column number of the 

end c* the Kth boundary sequence in the window. 

MINL MINL is a dummy variable which changes and represents the 

smaller of the window or picture boundary sequence minus one. 

CALL TIMNOW ( ) Subroutine for calling internal clock to time computation of 

NC’OR(U). 

The first step of the program is to read in the sections of the binary boundary 
maps to be considered, skipping all of the nonboundary or zero elements, and storing the 
start scan, start column, and column length of the sequences in the window and picture. 
This information is stored consistent with the definition of IW(K), JVV(K), NSQW(K), 
IP(L), JP(L), NSQP(L), KMAX, LMAX1, and LMAX, and results in a considerable 
compression of the data. The calculation of the above variables is accomplished in the 
DO 20 nested do-loop for the window and the DO 40 nested do-loop for the picture. 
The remaining part of the program contains two major calculation segments, the DO 100 
nested do-loop and the DO 200 nested do-loop. Because of their similarity it will suffice 
to explain only one of these do-loops and then point out their differences. The DO 100 
segment is concerned with the entire window and the top half of the picture; i.e., 
I = 1,2, ..., NROW. The two outer do-loops, DO 100 and DO 80, determine the 1 and J 
shifts that cause the individual boundaries of the K window sequences to coincide with 
the individual boundaries of the L picture sequences. No negative shifts are permitted; 
i.e., KKNCOL + l and KJ<NCOL+l. Within the DO 100 and DO 80 loops, 
there are 3 cases to be considered. The first case is the DO 56 loop and it occurs when 
the minimum length, MINL, of either the window or picture sequence is one boundary 
element. This occurrence is illustrated in Figure 3 and can result in an addition of at 
most 1 to NCOR(I.J) for each value of J. The case which probably occurs most often is 
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Figure 3. Subcase for DO 56 loop. 


the DO 68 loop, and a typical example is shown in Figure 4. The first time through the 
DO 68 loop NAD - 1 is added to NCOR(I,J) and NAD continues to increase by one until 
it is greater than M1NL, the minimum length of either sequence minus one. Once NAD is 
one greater than MINL it remains that value until J is greater than JM, and then NAD 
decreases by one each time until the two sequences no longer overlap in the column 
direction. Thus, both the initial and final values of NAD added to NCOR(l,J) are one, 
while the maximum value of NAD is the number of boundary elements in the smaller of 
the two sequences. Table 1 is a convenient way of viewing the number of occurrences of 
each column shift due to the presence of both sequences. The top two rows list the 
possible values of J or the column shift, while the next three rows list the possible values 
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Figure 4. Subcase for DO 68 loop, JL > 0. 
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TABLE I. NUMBER OF OCCURRENCES OF EACH COLUMN SHIFT 

FOR FIGURE 4 



JL 



JM 


JU 


3 

4 

5 

6 

7 

8 

JP(L) 

JL 

JL+1 

JL+2 

JL+3 



JP(L) + 1 


JL+1 

JL+2 
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JL+4 


JP(L) + NSQP(L) 



JL+2 

J'i-3 

JL+4 

JL+5 

NAD 

1 
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3 
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of J for each boundary element in the picture. The last row is the corresponding value 
for NAD for each J or the total number of occurrences of each column shift. An 
additional subcase for the DO 68 loop can occur when JL is negative, JP(L)> JW(K), 
and the length of the window sequence is equal to or less than the picture sequence. An 
example of this subcase is shown in Figure 5 which corresponds to Table 2. 

The final two subcases occur in connection with the DO 78 loop. Both subcases 
have the minimum column shift, JL, which is initially negative. The first subcase has 
JP(L) less than or equal to JW(K) and is illustrated in Figure 6 which corresponds to 
Table 3. The second subcase is illustrated in Figure 7 which corresponds to Table 4. 
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Figure 5. Subcase for DO 68 loop, JL < 0, NSPAN < NU, JP(L) > JW(K). 
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TABLE 2. NUMBER OF OCCURRENCES OF EACH COLUMN SHIFT 

FOR FIGURE 5 
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JM 



JU 
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Figure 6. Subcase for DO 78 loop, JL < 0, JP(L) < JW(K). 
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TABLE 3. NUMBER OF OCCURRENCES OF EACH COLUMN SHIFT 

FOR FIGURE 6 



JL 


JM 
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Figure 7. Subcase for DO 78 loop, JL < 0, JP(L) > JW(K), NSPAN > NU. 
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TABLE 4. NUMBER OF OCCURRENCES OF EACH COLUMN SHIFT 

FOR FIGURE 7 



JL 



JM 

JU 


1 
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4 

5 
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The second segment of the program or the DO 200 loop is identical to the first 
segment of the program, the DO 100 loop, except for the statement numbers. The main 
difference is that the DO 200 loop, while considering th,e entire window, only considers 
the bottom half of the picture. Thus, the DO 100 loop considers less and less of the first 
half of the picture sequences for all of the window sequences, ■ nd the DO 200 loop 
considers less and less of the window for the last half of all of the picture sequences. In 
the DO 100 loop, the scan shift I starts at the maximum available value and decreases to 
1=1, while in the DO 200 loop I starts at the minimum available value and increases to a 
maximum of I = NROW + 1 . A scan shift of I = 1 and a column shift of J = 1 indicate 
that the window and picture are already registered. 


IV. REGISTRATION SCHEME 


In applying the registration scheme it is assumed that the data can be misaligned 
by translation and/or rotation, but that no distortion exists between the two digital 
images. If significant distortions do exist over the entire images, then it may become 
necessary to register subportions of the digital images as separate data sets. 

It is also assumed that both digital images are of the same scale. If the images are 
not of the same scale, then it is possible to scale down one image by locating several 
corresponding boundary elements on both images, choose a corresponding reference 
boundary element on each image, and compute for each image the average distance of all 
the other boundary elements from the reference boundary element. The ratio of the two 
average distances gives a scaling factor. The raw data image with the larger average 
distance can then be scaled to the same size as the other image, and a new boundary map 
can be computed. 
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To register two boundary maps, choose one map as the reference map or picture 
and choose subportions of the other map to be used as windows as illustrated in Figure 
8. Figure 9 shows a portion of the two maps overlayed and the minimum and maximum 
scan and column shifts for the window. For the registration to work properly, the final 


IMAGE TO BE REFERENCE IMAGE 

REGISTERED OR PICTURE 



Figure 8. Window selection. 


(l-NROW+1, 



Figure 9. Overlay of digital images. 





two ; ssumptions are that the portions of the picture corresponding to the window are 
within the allowable range of the window and that the amount of rotation present in the 
picturr or the window can be neglected as far as the computation of NCOR(l,J) is 
concerned and can be treated as a local translation. If the latter assumption is not valid, 
it may be necessary to' repeat the registration procedure again after the data have been 
trams', ted and/or rotated based upon the first registration effort. In the event that the 
first assumption is correct, it should be possible to overcome the second assumption by 
repe >ied application of the registration procedure to the data. 

The mathematics of the translation-rotation required to register the two images is 
} .ivei below. Let w xj, w yj be the scan and column coordinates respectively of the upper 

left comer on the ith window in the image to be registered, and let lj and Jj correspond 

to the maximum value of NCOR(l.J) for the ith window. For any one of the windows, 
for < '.ample the jth, let the coordinates 


P x j = 


w x j 




and 


( 4 ) 


P y j = 


wVj 


+ J: - 1 


'>n the picture be designated as the center of rotation. Thus, if the entire window image 
.1 shifted by the amounts lj-1 and Jj-I, the coordinates ( p Xj, p yj) and ( w xj, w yj) will 

coincide and the rest of the image can be rotated about these coordinates to complete 
the registration, rwided a rotation exists. If no rotation exists, the images are registered 
and the maximum value of NCOR(l.J) for the rest of the window occurs at 1 = J= I. 
Assuming that rotation does exist, it will be necessary to compute the following 
quantities. 


1 

A* = S7, S VH ' pKj><w*i - w*j> 
r*3 


P y w y " N-| 



pyj) <wVi ~ wVj) 


16 ) 


13 


* 

k 





pXj) ( w Yj " ^j) 


( 7 ) 


pW 


_1_ 

N-l 


N 

S (p*i - 

B*3 



1 N 

E/w x i “ w^j) (py» - pyp 

B*3 


and 


( 8 ) 


p^y ” w*iP 
pV + 


( 9 ) 


where N is the total number of windows. To fetch the data from the image to be 
registered, it is first necessary to solve for 0, and for each picture coordinate ( p x, p y), it 
is necessary to solve the equations, 


p x cos0 + p y sin0 - p Xj cos0 - p yj sin0 + p Xj = w x 


and 


(10) 


- p x sin0 + p y cos0 - p xj sin0 - p yj cos0 + p yj = w y 


for the appropriate unregistered image coordinates. In general, the coordinates ( w x, w y) 

will not be integers, but in most cases it is probably most appropriate and less time 
consuming to round off to the nearest integer coordinate rather than try to interpolate 
the value of the data between the original data points. 

The registration routine is currently being programmed and debugged by MSFC’s 
Computation Laboratory, Engineering Computation- Division (S&E-COMP-RRV). The 
routine has not yet been modified to accept the correlation scheme, but instead accepts a 
manual input of corresponding boundary elements from two different images. 
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V. TEST CASES, RESULTS, AND SUMMARY 


The binary correlation routine was compared to a correlation and Fast Fourier 
Transform (FFT) correlation routine by obtaining running times from an internal clock 
on the IBM-7044 computer. A listing of the correlation routine and a piogram for using 
the FFT correlation routine are presented in Appendices B and C. The Fast Fourier 
Transform used, FFT7, is not presented in this report since it has been evaluated and 
listed in Reference 4 along with nine other FFT routines. According to Reference 4, 
FFT7 was the most versatile and ranked third in speed to the fastest routine, FFT3, 
which was twice as fast but worked only on real data, performed only a (+i) transform, 
and produced only a half transform. The second fastest transform was FFT8. which was 
only 25 milliseconds faster than FFT7. 

The data set on which the programs were run is shown in Figure 10. In the data 
set, only the locations of the boundaries are printed out and represented as I’s, and the 
window and the picture both start in the upper left comer. The size of the data array 
used for the window is the First NCOL columns by NROW scans, and the picture is twice 
as large. For the Fast Fourier Transform the window has to be the same size as the 
picture, and consequently zeros had to be added to the window array. The window and 
picture data arrays were arranged such that they were already in registration or such that 
the maximum value of NCOR(I,J) occurred at I = J = 1 . 

The output of the binary correlation routine for a 32 by 32 window array is 
shown in Figure 11, and the maximum number of NCOR(l,J) is the total number of 
boundary elements in the window. Figure 12 is a graph of correlation lag array size 
versus running time for the binary correlation. Fast Fourier Transform correlation, and 
correlation routines. The actual running times for each program versus correlation array 
size are listed in Table 5. Storage requirements for the FFT correlation routine were such 
that it was not possible to compute 41 by 41 correlation lag points or larger, and the 
irregular running times as shown in Figure 12 are a result of the FFT routine being more 
efficient for some array sizes than others. The correlation routine was able to handle 57 
by 57 correlation lag points, but the running time was already 5 minutes for 41 by 41 
correlation lag points. Since it appears that running times for the correlation routine can 
easily be extrapolated from Figure 12 and Table 5. no larger size arrays were run. 

According to Reference 1, the Sequential Similarity Detection Algorithm (SSDA) 
mentioned therein is approximately 50 times faster than the FFT correlation method. 
The factor of 50 is estimated from time ratios relating arithmetic operations on the IBM 
360/65 used in computing the numerator of the correlation coefficient. By relating all 
arithmetic operations to integer adds for a 32 by 32 window and a 256 by 256 picture, 
the equivalent integer adds for the correlation, FFT correlation, and SSDA respectively 
were 2.57 X 10 8 and 2.2 X 10 6 , which indicate a factor of 100 times faster than the 
FFT correlation. To compare the binary correlation with the SSDA it was necessary to 
put a counter in the DO loops containing NCOR(l,J) and to use the formula of 
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Figure 1 2. Correlation lag array size versus running time. 

TABLE 5. RUNNING TIME VERSUS CORRELATION LAG ARRAY SIZE 


Lag 

Array Size 


25 by 

25 

29 by 

29 

33 by 

33 

37 by 

37 

41 by 

41 

45 by 

45 

49 by 

49 

53 by 

53 

57 by 

57 

81 by 

81 

101 by 

101 


Running Time in 60' 1 seconds/Times Slower Than 
Binary Correlation 

Binary Correlation 

FFT Correlation 

Correlation 

58 

585/10.09 

2 450/42.24 

124 

1267/10.22 

4 465/36.01 

230 

983/4.27 

7 532/32.75 

381 

1778/4.67 

11 962/31.4 

559 
850 
1 167 

1 710 

2 377 
8 371 

19 761 


18 097/32.37 


i 











Reference 1 for computing the equivalent number of integer adds for a 32 by 3 2 window 
and 64 by 64 picture. The formula presented in Reference 1 is 4(1 + 10-y/\4/32XL * M+1) J 
equivalent integer adds, where L and M are the picture and window length or width, 
respectively. For the above mentioned picture and window the number of equivalent 
integer add* for the SSDA is 47 f 916 while the number of integer adds for the particular 
data set shown in Figure 10 is 20,804. Thus, it would appear that the binary correlation 
could possibly be 2.3 times faster than the SSDA. However, it is also possible that 
computing equivalent integer adds for various types of programs does not give 
the complete story, since it is necessary in all the computer programs to execute 
additional computer statements other than updating the correlation. Otherwise, the 
binary correlation would be approximately 200 times faster than the FFT correlation. 
Table S indicates that it is not. It is realized that one does not obtain “something for 
nothing,” and the processing time needed to produce a binary boundary map could 
negate the speed of the binary correlation. However, two points in addition to those 
mentioned in Section I are worth considering. According to Reference 1 it may become 
necessary to resort to the use of boundary maps when the channels under consideration 
are widely separated in spectral wavelength. This is certainly evidenced when one is 
working with near infrared or thermal imagery which tends to be quite noisy. Secondly, 
Figure 13 suggests a way to optimize the production of a boundary map. 

Figure 13 illustrates the number of variables needed to define various window 
sizes for the correlation and binary correlation methods using the data ret illustrated in 
Figure 10. Also included in Figure 13 is the number of vari? K s needed to define the 

location of the boundary elements only in the data set, which is twice the number of 

boundary points. For the correlation it is necessary to store the entire window array, 

while for the binary correlation it is necessary only to store 3 X KMAX variables, the 

start scan, the start column, and the sequence length for each boundary sequence. Thus, 
Figure 13 indicates the amount of data compression that is possible. Table 6 lists the 
compression factor for various window sizes, which is the total number of data points 
divided by 3 X KMAX. 

The production of a boundary map could be optimized by recognizing that three 
variables are needed to define each boundary sequence. Also, from the logic used in 
calculating IW(K), JW(K), NSQW(K), IP(L), JP(L), and NSQP(L), it is only necessary to 
determine boundaries within a scan and not necessary to determine boundaries that occur 
across two adjacent scan lines. Thus, the boundary program described in Section 11 could 
be optimized by working with one scan of data at a time instead of two and by 
eliminating all boundary sequences that contain less than three boundary elements. This 
would also optimize the compression by assuring that 3 X KMAX is always equal to or 
less than the number of boundary elements. 

It is foreseen that in some cases it may become necessary to subtract mean values 
from NCOR(M) and normalize to produce a proper registration of two data sets. In this 
respect it is interesting to examine the correlation coefficient in terms of rewards and 
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TOTAL NUMBER OF VARIABLES 
NEEDED TO DEFINE DATA SET 


FOR CORRELATION A 


2 X NUMBER OF x 

BOUNDARY POINTS • 3 X KM AX 


TOTAL NUMBER OF DATA POINTS 

Figure 13. Number of variables needed to define window array. 


penalties used in the binary correlation method. The binary correlation method presented 
in Section III only rewards for boundary elements that match and does not penalize for 
mismatches of boundary and nonboundary elements. By subtracting mean values and 
normalizing, it is possible to obtain rewards, for matching boundary elements and 
matching nonboundary elements, and penalties for mismatches of boundary and 
nonboundary elements. Before examining the correlation coefficient, it will be necessary 
to define a few terms. 




TABLE 6. COMPRESSION FACTOR FOR WINDOW ARRAY 


Window 
Array Size 

Number of 
Boundary Elements 

KM AX 

Compression 

Factor 

24 by 24 

61 

33 

5.81 

28 by 28 

90 

48 

5.44 

32 by 32 

122 

70 

4.88 

36 by 36 

158 

95 

4.55 

40 by 40 

183 

115 

4.64 

44 by 44 

237 

151 

4.27 

48 by 48 

296 

178 

4.31 

52 by 52 

368 

229 

3.94 

56 by 56 

450 

281 

3.72 

64 by 64 


355 

3.85 

72 by 72 


451 

3.83 

80 by 80 

903 

552 

3.86 

88 by 88 


656 

3.93 

96 by 96 


789 

3.89 

lOOby 100 

1431 

847 

3.94 


NWB The number of boundary elements in the window. 

NPB(M) The number of boundary elements in the picture that corresponds to 

the area covered by die window. NPB(U) chances as the window is 
moved to different 1 and J locations. 
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W The mean value of the data in the window, W*NWB/l (NROWXNCOL)] . 

P(l,J) The mean value of the data in the picture that corresponds to the area 

covered by the window, P(I,J)*NPB(I,J)/[(NROWXNCOL)] . 

The numerator of the correlation coefficient, C(1,J) can then be written as 
C(1,J) * K, • NCOR(l,J) + Kj * [NPB(U) - NCOR(U)] 

+ K 3 • [NWB - NCOR(U)l + K 4 * [NROW • NCOL - NPB(I,J) 
+ NCOR(U) - NWB) , (.1) 

where 

K, = (1 - WJU - P(U)1 . Kj * -W[l - P(I,J)], Kj = -H-WIP(U) 


K 4 = W • P(I,J) 


Since 0<W, P(I,J)<1, then K ( ,K 4 >0 and K 2 ,K 3 <0, and K, and K 4 are always 
greater than K 2 and K 3 in absolute magnitude. If the constants, K, are thought of in 
terms of rewards and penalties, then the first term of equation (1 1) is a reward times the 
number of boundary elements that nutch in the window and picture. The second and 
third terms, since they are negative, penalize the mismatch of boundary and nonboundary 
elements in the picture and window, respectively. The last term is a reward for 
nonboundary elements that match in the picture and window. It is inteiesting to note 
that for W+P(I,J)<1, then K, > K 4 > K 2 ,K 3 ; that for W+P(U)=1. tlvn 
K| * K 4 > K 2 ,K 3 ; and that for W+P(1,J) > 1 , then 1^ > K, > K, . This indicates that 
the binary correlation procedure gives a larger reward for matching up whatever occurs 
the least in the window or picture, boundaries or nonboundaries. This is also what a 
human observer would tend to do. To determine a normalization factor for equation 
(11), consider a perfect match for the window and picture; i.e., NWB*NPB(U )*NCOR(U ) 
and W»P(I,J) < 1. Equation (11) becomes 
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C(1J) * (1-W) 2 • NWB + W 2 • (NROW • NCOL - NWB) 


- 11 -P(U)] 2 NPBOJ) + IP(U)1 2 

INROW- NCOL-NPB(U)l (12) 


Thus, for a nonperfect match, the normalization factor should be some combination of 
the two expressions in equation (12). Since the two expressions in equation (12) are the 
variances of the window and picture respectively, the correct normalization factor would 
be the square root of the products of the two expressions to produce a correlation 
coefficient. 

It does not appear that computing the correlation coefficient would add a 
significant amount of running time to the binary correlation routine. This is 
because the mean value and variance of the window need only be calculated once and 
KMAX 

NWB=KMAX + Yj NSQW(K). The calculation that would require additional time and 

K=1 

storage is determining NPB(I,J) for the mean value and variance of the picture for each 
shift of the window. In most cases, however, it is anticipated that it will be necessary to 
compute only NCOR(U)- 
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APPENDIX A. BINARY CORRELATION LISTING 


5245. CAMPBELL .140145,30. FORTRAN SOURCE LIST 

ISM SOURCE STATEMENT 

0 UbFTC MAIM DECK 

1 DIMENSION NNI200) ,FMT( 12) *NCOft 1 101,101). Ml 13241 »J4t 1324) 

2 DIMENSION NSQMI 1324) . I PI 3300 ) » JPI 33001 .NSQPt 3300 1 

3 READ 330. NROW.NCOL 


4 

300 

FORMATI2I51 

7 


PRINT 502.NR0H.NC0L 

10 

302 

FORMAT ( IX . 5HNR0W* *14.1 X.5HNC0L 

11 


READ! 5.301 MFMTILI.L-l.l2l 

16 

301 

FORMAT! 12A6) 

17 


MROWl-NftOW+1 

20 


MC0L1 *NCOL«-l 

21 


MC0RINR0M1 . NCOi.il «0 

22 


K«3 

23 


00 20 1»1.NR0W 

24 


IN»-1 

25 


READ! 5«FMT)(NW( 11).11<1.NC0L) 

32 


00 10 J«1.NC0L 

33 


MCORI I.JKO 

34 


IFINWt JUcQ.OIGO TO 5 

37 


IF(IN.GT.C) GO TO j 

42 


K-K+l 

43 


IWIK1-I-1 

44 


JMIKI *J— 1 

*5 


NSQWIKKO 

46 


IN®1 

47 


GO TO 10 

50 

3 

MSQW<KI<NSC0<KK1 

51 


GO TO 10 

52 

5 

IN«-1 

53 

10 

CONTINUE 

55 

20 

CONTINUE 

57 


KMAX«K 

60 


RE AOI 5,301 ) IFMTILl .L*l»12 1 

65 


MR0W2»MR0W«’NR0W 

66 


MC0L2*NC0L»NC0L 

67 


l»C 

70 


DO 40 I«l,NR0W2 

71 


In*- l 

72 


RE AOI 5»FMT MNWIIll.1 1»1 *NC0L2) 

77 


DO 30 J*1 ,NC0L2 

100 


IFINWI J).fc«.0)GJ TO 25 

iu3 


IF(IN.GT.C) GO TO 23 

106 


L»L*1 

107 


NSQPILKO 

110 


IPILKI 

111 


JPILKJ 

112 


IN»1 

113 


IF 1 I . GT.NROMl I GO TU 30 

116 


LMAXl «L 

117 


30 TO 3C 

120 

ii 

NSQPIL KMSuPILKl 

121 


3U TO 10 

122 

25 

|N«-1 

123 

»Ki 

C'J'lTINUr. 

125 

>0 

CONtlNU- 
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5**5. CAMPBELL , 1*01*5, 30, FORTRAN SOURCE LIST MAIN 


1SN 


SOURCE STATEMENT 



127 


LMAX-L 



130 


XSTRT-l 



131 


I SUN-0 



132 


CALL T1MN0WUT1) 



133 


00 100 L-l.LMAXl 



13* 


NJ»JP(U*NS9P(L» 



135 


00 80 K-l.KMAX 



13b 


lFIl«|K).GE.IP(L)>GO to 

100 


1*1 


I-IP(U-lM<K) 



1*2 


IF< I.Nt.lJGO TO 65 



1*5 


IHNU.GT.JW(K) )G0 TO 66 



W 


GO TO 100 



151 

*5 

IFINU.LE.JWtKDGO TO 80 



15* 

*6 

JL-JP ( L)**<JW(K )— NSQMt K) 



155 


1F1JL.GT.NC0L1IG0 TO 80 



160 


JU-NU-JWIK) 



161 


IFUL.LT.1IG0 TO 70 



16* 


IFtNSUPtU.LE.NSUWUnGO 

TO 

53 

l«»7 


MINL-NSQW(K) 



17U 


GO TO 56 



171 

53 

HINL-NSQPIU 



172 

5* 

IF (MINL.GT.O)GO TO 60 



175 


IFIJU.LE.NC0L1IG0 TO 55 



200 


JU-NCOLi 



20 1 

55 

00 56 J-JC.JU 



202 


I SUM* I SUM* 1 



203 

56 

NCORI I .JI-NCORII , Jl+l 



*05 


GO TO 30 



206 

oO 

JM»JU-N!NL 



207 


NAO-O 



210 


IFtJU.LE.NCOLDGO TO 66 



213 


JU-NCOLI 



21* 

66 

00 68 J-JL.JU 



215 


I SUM- 1 SUM+l 



216 


1FIJ.GT. JMJGO TO 66 



221 


IFNAO.GT. HINDOO TO 68 



22* 


NAD-NAO+l 



225 


GO TO 66 



22b 

66 

NAO-NAU-1 



227 

66 

NCORI I ,J)*NCORI I , JI-fNAO 



<.31 


SO TO SO 



232 

70 

NSPAN-JWlKl+NSlttf (lO + l 



233 


IF(JP(U.GT.JW(K))GO TO 

71 


236 


IFINSPAN.GE.NUIGO TO 69 



2*1 


NAO-rtSQMIK !♦! 



2*2 


GO TO 75 



2*3 

6V 

NAl.-JU 



2** 


SO TO !■> 



<.*5 

7: 

IFINSPAN.GE.NUISU TO 7i 



250 


NAO-NSPAN-JPIU 



251 


IFINSOP(L) .GT.NStiwUO (GO 

Til 

76 

-56 


MINL-NSQP(L) 



<;5i 


GO TO 76 



256 

76 

MtNL-NSUW(K) 



257 

76 

JM-JU-MINL 




25 


52*6. CAMPdtLL ,140145,00, FORTRAN SOURCE Lf'.T MAIN 


tSN 


SOURCE STATEMENT 

1 00 


JL-1 

lb i 


GO TO 6 * 

coc 

7 J 

NAU*N SQP 1 L !♦ 1 

263 

75 

JM*JU-NAU*2 

264 


00 78 J*1 » JU 

<65 


ISUM-ISUM+l 

<. 6 o 


IF ( J.LT. JM)GO TO 76 

271 


NAfJ*NAO-i 

272 

78 

NCOR t 1«J)*NC0R(1 , J ) +NAO 

274 

80 

CONTINUE 

276 

100 

CONTINUE 

300 


LMAX1*LMAX1*1 

301 


00 200 L-LMAAi.LMAX 

302 


NU«JP IL ) ♦NSQPI LI 

303 

52 

OU 160 K*KSTRT,KMAX 

30* 


I-1PILI-IMIK) 

305 


IF ( 1 . GT.NR0W1 ) GO TO 85 

310 


IF,NU.LE.JM(KI)GO TO 180 

313 


JL*JPIU-JW(K)-NSQNllO 

314 


IF( JL .6T.NC0L1 )G0 TO 180 

31? 


JU*NU“JW(K ) 

.>20 


IFIJL.LT. 1 160 TO 701 

323 


IFINSQPIU.LE.NSUMIKHGO TO 

526 


MINL*NSa*(IO 

327 


GO TO 5*1 

330 

531 

MINL-NSQPILI 

351 

5*1 

IFIMINL.GT.OI60 TO 601 

3 3* 


IFIJU.LS.NCQLUGO TO 552 

317 


JU-NCOLl 

540 

552 

DO 561 J*JL* JU 

3*1 


ISGM*ISJM +1 

3*2 


NCORI I,J)*NC0R(I,J>+1 

34* 


GO TU 180 

345 

bOi 

JM-JU-MINL 

3*0 


N A0*0 

3*/ 


IFIJU.Lrf.NCJLDGO TO 6*1 

352 


JU*NCOLi 

353 

6*1 

00 661 J*JL, JU 

35* 


1 SUM* 1 SUM* 1 

->55 


IFU.GT.JMIGU TO 661 

jou 


IFIMAO.GT.MINUGO TO 661 

363 


N Al)*NAU+ 1 

36* 


GO TO \«l 

365 

6 bi 

NAO-NAO-l 

36 b 

68 1 

NCGRI l , J ) =NCOR ( I , J ) ♦NAD 

.'7<: 


GO TO lbC 

3U 

7v. i 

NSPAN*JW(K)+*<SJW(K) + 1 

372 


IFIJPID.jT.JwlKIIGU TO 711 

375 


IFINSPAM.GE.NJIGO TU 691 

** 0 


NAt>«NSCb<* > ♦ l 

*01 


GO TO 751 

2 

6 *» 

'4Alt*JU 

*>: 3 


GO TO 751 

*0* 

711 

IFINSPaN.Gc.NUI^U TO 731 

*07 


NA0*i'iSP4N- Ji’I L 1 
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5245 • CAMPBELL ,140145.00. FORTRAN SOURCE LIST MAIN 


ISN 


SOURCE STATEMENT 

410 


IFINSQPiL) •GT.NSQM(K) IGO TO 741 

413 


MINL-NStiPiL) 

414 


80 TO 761 

413 

741 

MINL-NSQNIK) 

414 

741 

JM-JU-MINL 

41? 

"•* * 

JL-I 

- 420 


80 TO 641 

' 421 

731 

NAU«NS0PIL)*1 

422 

751 

JM-JU-NAD+2 

423 


00 781 J-l.JU 

424 


1 SUM- 1 SUM* 1 

425 


IFIJ.LT.JMIGO TO 781 

430 


NA0*NAD-l 

431 

781 

NCOR l I . J » -NC0R I 1 . J » *NAD 

433 

180 

CONTINUE 

455 


80 TO 200 

434 

65 XSTRT-K*! 

•'^43? 

*- .. 

60 TO 52 

. " * 40 

*2 

© 

o 

CONTtNUE 

442 


CAU TIMN0UIIT2) 

443 


PRINT 5555, 1 SUM 

444 

5555 

FORMAT 1 5X.7HI SUM - ,171 

445 


4TIME-IT2-IT1 

446 


PRINT 318 

447 

316 

FORMAT I5X,50H**************************************************) 

450 


PRINT 317, JTIME 

451 

317 

FORMAT! /////, 5X, 15he LAPSED TIME - , I5,llH/oO SECONDS,//////) 

452 


PRINT 1100,KHAX,LHAX 

453 

1100 

FORMAT! 5X.7HKMAX - , I5,5X,7HLMAX - ,15) 

454 


PRINT 318 

455 


PRINT 999 

456 

999 

FORMAT! 5X.15HRE5IS4 COMPLETE) 

437 


JSTRT-l 

460 


JST0P-20 

4ol 

405 

PRINT 403 

462 

40 j 

FORMA T !1H 1 , 5X , 18HC0RRE LAT ION MATRIX) 

463 


DO 401 U-1.NR0H1 

464 


PRINT 402 , INCORI I1,J),J»JSTRT*JST0P) 

*»71 

402 

FURMAT! IX, 20! I4,1X) ) 

472 

401 

CONTINUE 

474 


IFiJSTOP.20.NCOLl)80 TO 406 

477 


JSTRT-JSTQP+1 

300 


JST0P-JSTUP*20 

501 


IK JSTOP.LS.NCOLl )80 TO 405 

504 


jstdp-nluli 

505 


00 TO 405 

306 

406 

C'INT I NUE 

507 


STOP 

»:u 


cNU 
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APPENDIX B. CORRELATION ROUTINE LISTING 


52*5, CAMPBELL .1*01*5, CO. PORT It AH SOURCE LIST 

1SN SOURCE STATEMENT 

0 UBFTC COR OE'K 

1 SUBROUTINE «OR(NP,NM, NCORI 

C********?********************************************************************** 

t THIS SU&kOUT I N£ CONTAINS THE STANOARU CORRELATION CALCULATION * 

C******************************************************************* ♦**•*•**•*•* 

2 DIMENSION NMI5A.SAI .NPt 112.112 1 .NCORI 5T.ST I.PNTI III 

3 CONMON/NAME1/ NCOL 

* C3MNON/NAME2/ NROM 

5 CONMON/NAME*/ FNT 

6 READ (5,5011 (FNT(L> ,L*l, 12) 

15 501 F0RMATI12AA) 

1* 03 10 11*1 .NROM 

15 AEADIS.FMT) (NNdl.Kl), K1*1,NC0U 

22 10 CONTINUE 

2* NROM*NROM*NROM 

25 NCOL*NCQL*NCOL 

2* READ (5,501) (FNT (LI »L*l , 12 ) 

53 DO 20 12*1, NROM 

3* READIS.FNT) (NF(I2,K2), R2-1.NC0L) 

*1 20 CONTINUE 

*3 NR0M*NR0M/2 

A* NCOL-NCOl/2 

*5 CALL TIMNOU(ITl) 

*6 NAOMI •NROM*! 

AT NC JLl*NCOL*l 

50 00 1 1*1 ,NROMl 

51 00 2 J-l.NCQLl 

52 NCORI l,J)»0 

53 DO 3 X*1,NAQM 

5* !l«I*k-l 

53 DO A L*1,NC0L 

5A Jl-JtL-1 

57 NCOR( I ,J)*NCOR( 1 , JI+NM(K,L)*NP( 11, Jl) 

60 A CONTINUE 

62 3 CONTINUE 

6A 2 CONTINUE 

66 l CONTINUE 

70 CALL TIMNOMI IT2) 

71 JTIME-1T2-IT1 

72 FAINT 318 

73 316 FORMAt(5X,50H*******A**A****A***A*A*AA*AA******A*******AA*A***A) 

7* FAINT 317, JTIME 

75 317 FORMAT!/////, 5X.15HELAPSED TIME * ,I5,UH/60 SECONDS,////// 1 

76 FAINT 318 

77 PAINT AON 

100 9AN FORMAT ( 5X, I5HC0A COMPLETE) 

101 RETURN 

102 SNu 


28 


r 


APPENDIX C. FAST FOURIER TRANSFORM CORRELATION LISTING 

9209, CAMPBELL ,100109,00, FORTRAN S3URCI LIST 

ISM SOURCE STATE MINT 

0 MBFTC MAIN 

1 OINENSION XNMRt60»B0)»XNNI (00«60l«XNPftt60,60),XMPI(S0« 64) 

3 OINENSION |A(A4,*AI 

3 READ 1.NR0M.NCQL 

* I FORMATUISI 

T CALL tN(XNPM*XNPI,NAON«NCOL) 

10 CALL INtXNNR,XNNI,NRON,NCOLt 

U NTOT-NAONONCOL 

12 CALL TIMMMIITII 

IS CALL FFTT(XNPR,XNPI,NTOT,NROU,NCOL,1I 

.10 CALL FFTTIXNFR.XNF1 tNTOTtNROH.NTOT.il 

IS CALL PFT7IXNNA*XNMI «NTOT,MtON»NCOL,l> 

IA CALL FFTTIXNUfttXNMl ,NTOTkNROM*NTOT*ll 

17 XNTOTMTOT 

20 00 7 1*1,NA0N 

21 00 • J-l.NCOL 

22 StORE>9iM(l t J>*XNMR(l»JUXNFIU»2)*xmttll,JI 

23 MOt 1 1 ,JI»IXNP1I l « J>**NMAt ItJl-XNMH I « JI*XNPM( I ,J> l/XNTOT 

20 XNPR(1,J»*ST0AE/XNT0T 

29 8 CONTINUE 

27 7 CONTINUE 

91 CALL FFTTIXNPRtXNPI »NtOT,NROM,NCOL*~ll 

32 CALL FFT7<mPA,XNPt«NT0T«NR0M,NT0T,~l) 

33 CALL TINN0MIIT2I 

30 1TINE-IT2-IT1 

39 MINT 10 

30 FAINT IStITINE 

37 13 FORNATtUa.SXtEHlTINE • .I9.3M/60) 

00 FAINT 10 

01 10 

10000I 

02 CALL OUTtXNFR'NAOtttNCOLl 

03 STOF 

00 ENO 
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